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Abstract 


In this paper, a two-phase non-isothermal PEM fuel cell model based on the previously developed mixed-domain PEM fuel cell model with a 
consistent treatment of water transport in MEA has been established using the traditional two-fluid method. This two-phase multi-dimensional 
PEM fuel cell model could fully incorporate both the anode and cathode sides, properly account for the various water phases, including water 
vapor, water in the membrane phase, and liquid water, and truly enable numerical investigations of water and thermal management issues with the 
existence of condensation/evaporation interfaces in a PEM fuel cell. This two-phase model has been applied in this paper in a two-dimensional 
configuration to determine the appropriate condensation and evaporation rate coefficients and conduct extensive numerical studies concerning the 
effects of the inlet humidity condition and temperature variation on liquid water distribution with or without a condensation/evaporation interface. 


© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


Numerical modeling and simulation has been an important 
tool for facilitating PEM fuel cell design and optimization. Many 
multi-dimensional PEM fuel cell models have been developed 
in the past decade, and they can generally be categorized into 
two methods, namely the single-domain [1—9] and multi-domain 
[10-13] approaches. A mixed-domain model, which maintains 
a consistent treatment of water transport in the membrane- 
electrode assembly (MEA), has been recently developed [14] 
and further extended to investigate the effects of the fully coupled 
transport phenomena on fluid flows, species distributions, water 
transport processes, and detailed cell performances in PEM fuel 
cells [15]. 

Significant research efforts in the field of PEM fuel cell mod- 
eling and simulation have presently been focused on two-phase 
transport phenomena and the drastic effects on water man- 
agement and cell performances. The two-phase PEM fuel cell 
models in the open literature account for the liquid water trans- 


* Tel.: +86 571 87953166; fax: +86 571 87953167. 
E-mail address: menghua@zju.edu.cn. 


0378-7753/$ — see front matter © 2007 Elsevier B.V. All rights reserved. 
doi: 10.1016/j.jpowsour.2007.03.012 


port using either the multiphase mixture (MÊ) model [12,16-20] 
or the traditional two-fluid model [9,21—24]. A common weak- 
ness in these models, however, remains in its inability to predict 
significant liquid saturation values in the porous gas diffusion 
layer (GDL) and catalyst layer (CL) and consequently the flood- 
ing effects in PEM fuel cells. In the work of Meng and Wang 
[25], a liquid droplet coverage model at the GDL and gas chan- 
nel (GC) interface was proposed based on recent visualization 
experiments [26,27], and it was applied to successfully predict 
liquid water flooding dynamics. 

Although progress has been made in multi-dimensional mod- 
eling of two-phase transport phenomena and flooding effects 
in PEM fuel cells, significant improvements are still needed, 
particularly in the areas of fully integrating two-phase flows 
with heat transfer phenomena and properly accounting for the 
effects of low-humidity inlet conditions on liquid water trans- 
port and flooding dynamics, and under both circumstances 
could exist the condensation/evaporation interfaces. These are 
remaining crucial research areas to be properly handled before 
the intricate interplays of water and thermal managements in 
PEM fuel cells could be fully investigated. Recently, Ju et 
al. [28] and Luo et al. [29] presented a three-dimensional 
two-phase isothermal PEM fuel cell model for investigating 


H. Meng / Journal of Power Sources 168 (2007) 218-228 219 


Nomenclature 

Nomenclature 

a water activity or stoichiometry coefficient 

c molar concentration (mol m~?) 

Cp constant-pressure heat capacity (J(kg K)~!) 

D mass diffusivity (m? s7!) 

Dy water content diffusivity (mol(m s)7!) 

EW equivalent weight of the membrane (kg mol™!) 

F Faraday constant, 96,487C mol! 

hpc condensation/evaporation parameter 

hfg heat of vaporization (J kg7!) 

i current density vector (A m7?) 

j transfer current density (A m7?) 

k thermal conductivity (W(m K)— 1) 

ke condensation rate coefficient (s7!) 

ke evaporation rate coefficient (s7! Pa7!) 

K permeability (m7) 

M a symbol representing species 

na electro-osmotic drag coefficient 

p gas-phase pressure (Pa) 

Pe capillary pressure (Pa) 

Ru universal gas constant (J(mol K)"!) 

s liquid saturation 

S source term 

T temperature (K) 

u gas-phase velocity (m s7!) 

Uo open-circuit potential (V) 

Veell cell voltage (V) 

W molecular weight (kg mol~!) 

Greek symbols 

X mole fraction 

E porosity 

Em fraction of the membrane phase in the catalyst 
layer 

@ phase potential (V) 

n over-potential (V) 

K proton conductivity (S m7!) 

À water content 

u viscosity (kg(m s)7!) 

Ac contact angle 

p gaseous density (kg m7?) 

o electronic conductivity (S m~ hy or surface tension 
(Nm!) 

T viscous stress tensor 

Supercripts 

cl catalyst layer 

eff effective value 

l liquid phase 

sat saturation value 

v vapor phase 


Subscripts 

cl catalyst layer 

e electrolyte or energy 
gaseous phase 

i species 

int interfacial value 

l liquid 

m membrane 

s electron 

sat saturation value 

w water 


the condensation/evaporation interfaces under low-humidity 
inlet conditions, and they successfully predicted the impor- 
tant feature of dry-wet-dry transition in PEM fuel cells. 
Non-isothermal two-phase calculations have also been pro- 
vided in Ju [30]. All these simulations are based on the M? 
model. 

In this paper, the pseudo single-phase mixed-domain PEM 
fuel cell model with a consistent treatment of water transport in 
MEA [14,15] is further extended to consider two-phase flows 
and liquid water transport based on the traditional two-fluid 
model. The important features of the present multi-dimensional 
two-phase mixed-domain PEM fuel cell model are mainly in 
the following areas: (1) consistent treatment of water transport 
in MEA with the presence of liquid water and thus appropri- 
ately including both the anode and cathode sides in the model, 
(2) proper account of the various water phases and their trans- 
port processes in a PEM fuel cell, including water vapor, liquid 
water, and water in the membrane phase or the dissolved water 
phase, (3) incorporating the liquid droplet coverage model at the 
GDL and gas channel interface [25] and thus capable of exam- 
ining flooding dynamics and its effects on cell performances, 
and (4) fully coupling heat transfer phenomena with two-phase 
flows and consequently truly enabling numerical investigations 
of water and thermal management issues and their intricate 
interactions with the existence of a condensation/evaporation 
interface. 

This two-phase model is applied in this paper in a two- 
dimensional configuration to determine the key parameters of 
condensation and evaporation rate coefficients and also to enable 
extensive parametric studies, focusing mainly on the effects of 
the inlet humidity condition and temperature variation on liquid 
water distribution with or without a condensation/evaporation 
interface. 


2. Theoretical formulation 


The complete conservation equations of this multi- 
dimensional two-phase mixed-domain PEM fuel cell model are 
presented in this section. First, the conservation equations of 
mass, momentum, and species concentrations in the gaseous 
phase are established. 
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Mass conservation: 


V- (pu) =0 (1) 
Momentum conservation: 

— V -(puu) = —Vp+ V-1t4+58 2 

2 — s} (puu) P u (2) 
Species conservation: 

V : (ici) = V - (D**Ve;) + S; (3) 


In Eq. (2), the source term is added in the porous materials 
based on the Darcy’s law considering the liquid water effect, 


H 
4 
KK (4) 


Su = — 


The relative permeability is defined as [16] 


Ky =(1—s) (5) 


where the parameter, s, is the liquid saturation, defined as the 
ratio of the liquid volume to the pore volume [25]. 

In Eq. (3), based on the mixed-domain approach [14,15], the 
water vapor concentration is solved only in the gas channels, gas 
diffusion layers, and catalyst layers on both the anode and the 
cathode sides. Inside the membrane, a water content equation 
is solved, as presented later in this section. In the two cata- 
lyst layers, the dissolved water phase (water in the membrane 
phase) is assumed to be in thermodynamic phase equilibrium 
with water vapor, and its transport process is considered based on 
the “fictitious water concentration” treatment [14,31,32] using 
the following water diffusivity: 
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The water vapor and the dissolved water phase in the cat- 
alyst layer and GDL could interact with liquid water through 
condensation/evaporation processes. 

In the present two-phase model, water produced in the 
cathode catalyst layer is assumed to be in vapor phase as in 
[33] to be consistent with our previous pseudo single-phase 
method [14,15]. This is different from the assumption made in 
(21—24,34], but it should make no difference in the final results, 
i.e. the liquid saturation and temperature distributions, since the 
final states of water vapor and liquid water should approach 
their thermodynamic equilibrium conditions through conden- 
sation/evaporation processes in both methods, as discussed in 
detail later in this section and in Section 3. The source term in 
Eq. (3) can be expressed as 


(a) for species except water 


= —— 8 
nF a) 
(b) for water vapor 
ndz dw] 
Sy = -V (Si) - SF — sy (8b) 


In deriving Eq. (8) the electrochemical reactions are formu- 
lated in the following general form: 


y aiM; = ne 
i 


In Eq. (8b), the variable, Syı, is the volumetric condensa- 
tion/evaporation rate, which is defined as 


(9) 


RyT dr sat 
1 1.5 pel 1.5 u = Y 78a 
Dy, = eg? DYE + ex De a, (6) Sv = hpp =p) (10) 
Psat Ga 
PPE ee : where the parameter, Apc, is determined as [34 

Considering the liquid water effect, the effective gaseous P PEER [a 
species diffusion coefficients are further modified as at 

: p = keel = sity |, [PY = pl 

i 1.5 = 5 
Di" = Di — s) (7) ps 2RaT pY — pt 
Table 1 
Electrochemical and physical relationships 
Description Expression Unit 
CH: a CO: 
Transfer current density j= aji ( a) (gg SFe n) in anode side, j = aje! (<2) exp ( oe oe n) in cathode side Am? 
Over potential n = os — de in anode side, n = @s — de — Uo in cathode side V 
Open-circuit potential Uo = 1.23 — 0.9 x 107°(T — 298) V 
; - 1.0 fora < 14 
Electro-osmotic drag coefficient ng = { 1.5/8(à — 14) + 1.0, otherwise 
Water activity = oe 
Water saturation pressure logy) pt = —2.1794 + 0.02953(T — 273.15) — 9.1837 x 107$(T — 273.15)? + 1.4454 x 1077(T — 273.15) atm 
ees 3.1 - 1077a(e9:284 — 1) el-246/T] QQ <A <3 Szi 

Membrane water diffusivity DY = 417. 1078101 4 161e-%) el-2346/7 aea: mî s 
Water content diffusivity D, = & DP mol m7! s7! 
Proton conductivity k = (0.51391 — 0.326)exp [1268 (sas =d )] Sm! 


T. 
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The condensation and evaporation rate coefficients in Eq. 
(11), ke and ke, will be discussed and determined in the next 
section. The partial pressure of water vapor can be calculated as 


p` = CwRuT (12) 


An expression for calculating water saturation pressure is 
provided in Table 1. 

All the other relevant physicochemical relations and proper- 
ties can also be found in Table 1. 

In the present two-phase model, the effects of liquid water 
on the gaseous species transport are taken into account using 
Eq. (7) in the porous GDL and catalyst layer. In fact, inside the 
catalyst layer, with the presence of the membrane phase, liquid 
water distribution and its effects on species transport and the 
electrochemical reactions involve complex physical processes, 
which are not clearly understood. In the open literature, there 
were attempts to consider liquid coverage on the membrane 
phase and the ensuing extra resistance exerting on species trans- 
port [34,35], i.e. the agglomerate approach, but they generally 
applied many simplifications and thus caused uncertainty in the 
model accuracy. In the present numerical studies, we neglect 
this phenomenon and focus mainly on the effects of the inlet 
humidity and temperature variation on liquid water distribution 
with condensation/evaporation interfaces in the porous mate- 
rials. Furthermore, another reason for simplifying the catalyst 
layer model will be discussed in the next section when the 
numerical results are presented. 

In this two-phase model, liquid water transport is taken 
into account using the traditional two-fluid method to add the 
flexibility for studying the finite-rate condensation/evaporation 
processes. As discussed in the next section, we will ensure the 
states of water vapor and liquid water closely approach their 
thermodynamic equilibrium conditions by choosing appropri- 
ate condensation and evaporation rate coefficients in this paper. 
It should be noted that the two-fluid method is mathematically 
equivalent to the multiphase mixture M3 model [36]. 

Next, the conservation equations of liquid water are pre- 
sented. 

Liquid mass conservation: 


V. (pit) = Sy Ww (13) 


where an expression for the condensation/evaporation rate 
is presented in Eq. (10). In the porous materials, based on the 
Darcy’s law, the liquid water velocity is derived as 
> KaK 1 


u; = ———_Vp 


z 14) 


where the liquid pressure and the relative permeability are 
defined as [16] 


P =p- p (15) 


and 
Ka=s° (16) 


Eqs. (13)-(15) can be combined to produce a conservation 
equation for the liquid saturation, which is 


Kuk a KaK 
y. E ’ rys V. E r vp| =SuWw 
u Os H 


The capillary pressure, pc, can be further expressed as [16] 


p= (=) ETES (18) 


where J(s) is the well-known Leverett’s function, which takes 
the following form for a hydrophobic porous material with a 
contact angle larger than 90° [37]: 


J(s) = 1.417s — 2.120s* + 1.2635? 6, > 90° (19) 


Since numerical studies in this paper will be focused on two- 
phase transport phenomena in porous materials, liquid water 
transport in the gas channel and its interactions with the gaseous 
fluid-field are neglected. A mist flow model is used in Ref. [25]. 

Finally, the other conservation equations are derived as 

Proton transport: 


V - (kV pe) + Se =0 (20) 
Electron transport: 

V - (Vo) + Ss =0 (21) 
Energy: 

V - (PCpüT) + V+ (p1CptiT) = V-(RVT) + Sp (22) 


Water content conservation inside the membrane: 
V -(DyVA) + S, = 0 (23) 
In Eqs. (20) and (21), the source terms can be expressed as 
Se = —Ss = j (24) 


In the energy equation, Eq. (22), the effects of species dif- 
fusion on temperature distribution has been neglected, and the 
source term is given as [38] 


dUo i 
S =j (1 + ro) + a + Afg Ww Svi (25) 

where the second term inside the bracket is considered only 
on the cathode side. The parameter, cond, in Eq. (25) is either 
xt or of, depending on the location in the fuel cell. The last 
term in Eq. (25) is the heat of vaporization. 

In the water content equation inside the membrane, Eq. (23), 
the source term is 


gy. {ne (26) 
: F 


In addition, there are two boundary conditions at the inter- 
faces between the membrane and the two catalyst layers, which 
are needed to connect the water content and the water vapor 
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Cathode 


Anode 


Fig. 1. A two-dimensional cross section and the related boundaries. 


concentration equations [14,15]. The first boundary condition, 
assuming thermodynamic phase equilibrium at the interfaces, is 


[14,15] 
when s < 0 

à = 0.043 + 17.81a — 39.854? + 36.0a7 (27a) 
when s>0 

à = 1442.85 (27b) 


where the parameter à represents the water content on the 
membrane side of the interface while the parameters, a and s, 
the water activity and liquid saturation on the catalyst layer side 
of the interface, respectively. Eq. (27b) is proposed in Ref. [24]. 

The second boundary condition, ensuring water flux equality 
at the interfaces, is [14,15] 


naqi Leff nai 
(ov + i) = (-»% S VCw + ri) 


where the effective water diffusivity in the two catalyst layers 
is provided in Eqs. (6) and (7). 

The conservation equations, Eqs. (1)—(3), (17), and (20)-(23), 
constitute the present two-phase non-isothermal mixed-domain 
PEM fuel cell model. This model is applied for numerical sim- 
ulations in a two-dimensional configuration, as shown in Fig. 1, 
to determine the condensation/evaporation rate coefficients, and 
also to make extensive parametric studies. The computational 
domain includes five regions, namely the gas diffusion layers 
and catalyst layers on both the anode and cathode sides, and 
the membrane. In the present two-dimensional simulation, since 
fluid flows in gas channels are neglected, the gaseous velocity 
and pressure are not computed. This simplification is valid based 
on the theoretical and numerical analyses of the Peclet number 
in the porous materials in [15,39]. 

There are five boundary conditions to be specified, as shown 
in Fig. 1. Boundary 1 is at the GDL and gas channel interface 


(28) 


m cl 


on the anode side. The boundary conditions are defined as 


Ci = Cio (29a) 

S = Sint (29b) 

dhe Os 

dpe OPs _ 2 
ay Oe (220) 
T 

ar, (29d) 
Ox 


where the parameters, Cio and Sint, are specified values. In 
Eq. (29d), the heat flux is set as zero since the prior numeri- 
cal calculations indicate that the heat transfer rate through this 
interface is very small [38]. 

Boundary 2 is at the GDL and gas channel interface on the 
cathode side, the same type of boundary conditions as those in 
Eq. (29) are specified, and they are not repeated here. 

Boundary 3 is at the interface between the GDL and the 
current-collecting land on the anode side. The boundary con- 
ditions are defined as 


aC; 

Mi 0 (30a) 
Ox 

3 

29 (30b) 
Ox 

3 

dpe o (30c) 
Ox 

bs = 0 (30d) 
T= (30e) 


Boundary 4 is at the interface between the GDL and the 
current-collecting land on the cathode side. The same type of 
boundary conditions as those in Eq. (30) can be specified except 
for the one with the electronic phase potential, which becomes 


Ps = Veell (31) 


At boundary 5, the symmetric boundary conditions are used 
for all the variables. 


3. Result and discussion 


The present two-phase non-isothermal mixed-domain PEM 
fuel cell model has been implemented into a commercial CFD 
package, Fluent, through its user coding capabilities and applied 
herein for two-dimensional numerical simulations, as shown in 
Fig. 1. The geometric parameters of the fuel cell are listed in 
Table 2. 

The fuel cell is operated at 2 atm on both the anode and cath- 
ode sides. Hydrogen and water vapor is fed into the anode while 
air and water vapor into the cathode. For all the calculations car- 
ried out in this paper, the cell voltage, Vie, is fixed at 0.65 V, 
and the boundary temperature at 80°C. In order to investigate 
the effect of the inlet humidity conditions on liquid water dis- 
tributions, three cases have been designed, as listed in Table 3. 
Based on these operation conditions, the inlet species concen- 
trations can be easily determined and specified at boundaries 1 
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Table 2 
Cell geometric parameters 


Table 4 
Physicochemical parameters 


Fuel cell geometry [mm] 


Layer thickness 
Diffusion 0.3 
Catalyst 0.01 
Membrane 0.025 
Land width 0.5 
Channel width 1.0 


Computational cell numbers ~1600 


and 2. Therefore, the present two-dimensional numerical calcu- 
lations simulate a cross section perpendicular to the membrane 
at the inlet of a PEM fuel cell. The other relevant physical and 
transport parameters are provided in Table 4. Careful grid inde- 
pendence study has been conducted before extensive numerical 
investigations, and a total of 1600 computational cells have been 
used in the present calculations. 

Numerical studies are first conducted to determine an appro- 
priate condensation rate coefficient, ke, which is related to the 
surface to volume ratio in a porous material [33]. The calcula- 
tions are based on fully humidified inlet conditions, case 1 in 
Table 3, with a constant cell temperature of 80°C. This case is 
chosen because no evaporation process occurs and thus the evap- 
oration rate coefficient is irrelevant. Fig. 2 illustrates the liquid 
saturation distributions in cathode GDL and catalyst layer under 
three different condensation rate coefficients with a zero liquid 
saturation value specified at the boundaries 1 and 2. When the 
condensation rate coefficient increases from 500 to 2000 s7!, 
the maximum liquid saturation increases more than 20%, indi- 
cating the results are dependent on this parameter. Increasing the 
condensation rate from 2000 to 5000 s~! only brings around 6% 
increase of the maximum liquid saturation. Further increasing 
the condensation rate coefficient could hardly cause any increase 
of the maximum liquid saturation and variation of the liquid 
water distribution. The same conclusion could be drawn from 
Fig. 3, in which the water vapor concentrations under three corre- 
sponding conditions are presented. Fig. 3c shows clearly that the 
water vapor concentration approaches its thermodynamic equi- 
librium value at 15.9 mol m~, with only slight over-saturation. 
These results conclude that numerical results are independent 
of the condensation rate coefficient once it reaches a value 
of 5000s~!. Based on the present parametric studies and the 
detailed theoretical analyses in Ref. [33], a condensation rate 
coefficient at 5000 s~! is considered as an appropriate parameter 
and is thus chosen in the present two-phase PEM fuel cell mod- 
eling. Using this condensation rate coefficient, the calculated 


Anode volumetric exchange current density, ajo (A m~) 1.0E+9 
Cathode volumetric exchange current density, ajo (A m7?) 1.0E+4 
Reference hydrogen concentration, Cy, (mol m7?) 40 
Reference oxygen concentration, Co, (mol m~?) 40 
Anode transfer coefficients &a ==] 
Cathode transfer coefficient Qc = 
Faraday constant, F (C mol7!) 96487 
GDL porosity 0.6 
Porosity of catalyst layer 0.12 
Volume fraction of ionomer in catalyst layer 0.4 

GDL permeability (m?) 1.0E—12 
Catalyst layer permeability (m°) 1.0E—13 
Equivalent weight of ionomer (kg mol~!) 1.1 

Dry membrane density (kg m~?) 1980 
Electronic conductivity in current collector (S m7!) 20,000 
Effective electronic conductivity in GDL (S m7!) 5000 
Operation pressure (atm) 2 
Condensation rate coefficient (s7!) 5000 
Evaporation rate coefficient (s7! Pa~!) 1.0E—4 
Liquid water density (kg m~’) 1000 
Liquid water viscosity (Ns m~?) 3.5E—4 
Surface tension (N m~!) 6.25E—2 
Contact angle in GDL 110 
Contact angle in CL 95 
Thermal conductivity of current collector (W m7! K!) 20 
Thermal conductivity of the membrane (W m~! K7!) 0.5 

Heat of vaporization (J kg~!) 2.3E+6 


water vapor concentration and liquid water saturation approach 
the thermodynamic equilibrium condition. 

Next, parametric studies are conducted to determine an 
appropriate evaporation rate coefficient, ke. The calculations are 
based on case 3 in Table 3 under an isothermal condition at 80 °C. 
Liquid saturation distributions with condensation/evaporation 
interfaces under three different evaporation rate coefficients are 
illustrated in Fig. 4. Decreasing the evaporation rate coefficient 
by an order of magnitude, from 1 x 107° to 1 x 1074 (Pa s)—, 
brings no difference in the maximum liquid saturation value 
and the liquid water distribution around that region. However, 
decreasing the evaporation rate coefficient would result in a 
slightly larger liquid water distribution region. In the present 
two-phase model, an evaporation rate coefficient of 1 x 1074 
(Pas)~! is chosen as it makes the condensation and evaporation 
parts of the phase-change parameter in Eq. (11), 7pc, at the same 
order of magnitude. 

After the condensation/evaporation rate coefficients are 
determined, the two-phase model is then applied to study the 
effect of the liquid droplet coverage at the GDL and gas channel 


Table 3 

Inlet humidification temperature and relative humidity at 80°C (the fully humidified water vapor concentration is 15.9 mol m~?) 

Case number Anode Cathode 
Humidification Relative humidity Humidification Relative humidity 
temperature (°C) (%) temperature (°C) (%) 

Case 1 80 100 80 100 

Case 2 80 100 75 81 

Case 3 80 100 70 66 
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0.080 
0.075 
0.070 
0.065 
0.060 
0.055 
0.050 
0.045 
0.040 
0.035 
0.030 
0.025 
0.020 
0.015 
0.010 


Zz 


(a) (b) 


0.098 0.104 
0.092 0.097 
0.085 0.091 
0.079 0.084 
0.073 0.077 
0.067 0.070 
0.060 0.064 
0.054 0.057 
0.048 0.050 
0.041 0.044 
0.035 0.037 
0.029 0.030 
0.023 0.023 
0.016 0.017 
0.010 0.010 
Z 
X X 


(c) 


Fig. 2. Liquid saturation distributions at different condensation rate coefficients: (a) 500 s7}, (b) 2000 s7}, and (c) 5000 s~!. 


(GC) interface on liquid water distribution and cell performance 
through specifying different liquid saturation values at boundary 
2, a phenomenon also discussed in Ref. [25]. In Fig. 5, increasing 
the liquid saturation value at the GDL/GC interface on the cath- 
ode side, meaning more liquid droplet attached at the GDL/GC 
interface [25], the maximum liquid saturation in the GDL and 
catalyst layer on the cathode side increases drastically while 
the liquid saturation gradient decreases, a phenomena resulting 
from a much stronger capillary-driven flow in the porous mate- 
rials at a higher liquid saturation value. Similar results are also 
presented in Refs. [25,34]. The effect of the liquid droplet cov- 
erage at the GDL/GC interface on cell performance is depicted 
in Fig. 6, in which current density distributions along the lat- 
eral direction at three different liquid saturation boundary values 


19.30 
19.06 
18.81 
18.57 
18.33 
18.09 
17.84 
17.60 
17.36 
17.11 
16.87 
16.63 
16.39 
16.14 
15.90 


(a) (b) 


are presented. Increased liquid saturation values in the porous 
materials with increased liquid saturation boundary values on 
the cathode side reduce the cell performance. For example, 
increasing the parameter, Sint, from 0 to 0.35, decreases the cell 
performance by more than 10%. Furthermore, current density 
decreases more significantly under the land than under the gas 
channel, a result attributable to extra resistance for oxygen trans- 
port in the region. It should be emphasized that since this is a 
two-dimensional calculation without the contact resistance for 
the electronic transport, the cell performance obtained is thus 
higher than the practical value, but this will not affect the present 
numerical investigations. 

The effect of the inlet humidity conditions on the cathode 
side on liquid water distributions is clearly illustrated in Fig. 7, in 


17.90 17.30 
7.76 17.20 
7.61 17.10 
7.47 17.00 
7.33 16.90 
7.19 16.80 
7.04 16.70 
6.90 16.60 
6.76 16.50 
6.61 16.40 
6.47 16.30 
6.33 16.20 
6.19 16.10 
6.04 16.00 
5.90 15.90 

Z Z 
x X 


(c) 


Fig. 3. Distributions of water vapor concentration at different condensation rate coefficients: (a) 500 sL, (b) 2000 s~!, and (c) 5000 sv, 
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0.053 
0.050 
0.047 
0.044 
0.041 
0.038 
0.035 
0.032 
0.028 
0.025 
0.022 
0.019 
0.016 
0.013 
0.010 


(a) (b) 


0.053 0.053 
0.050 0.050 
0.047 0.047 
0.044 0.044 
0.041 0.041 
0.038 0.038 
0.035 0.035 
0.032 0.032 
0.028 0.028 
0.025 0.025 
0.022 0.022 
0.019 0.019 
0.016 0.016 
0.013 0.013 
0.010 0.010 
Z Z 
X X 


(c) 


Fig. 4. Liquid saturation distributions at different evaporation rate coefficients: (a) 1.0E—3 s7! Pa~!, (b) 5.0E—4 s7! Pa~!, and (c) 1.0E—4 s7! Pa™!. 


which results correspond to the three cases in Table 3 with a con- 
stant cell temperature of 80°C are shown. With a low humidity 
inlet condition at a humidification temperature of 70°C, liq- 
uid water only exists under the two land areas while the large 
area under the gas channel is at dry condition. The condensa- 
tion/evaporation interface is well captured and can be clearly 
seen in Fig. 7a. Increasing the inlet humidification temperature 
on the cathode side to 75°C results in drastic decrease of the 
dry region, and now the condensation/evaporation interface is 
completely under the gas channel, as shown in Fig. 7b. The dry 
region disappears under the fully humidified inlet condition as 
expected in Fig. 7c. Three-dimensional simulations predicting 
dry-wet-dry transitions in PEM fuel cells have recently been 
studied in Refs. [28,29]. 
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Fully integrating water transport, including possibly liquid 
water transport, with heat transfer phenomena is crucial for 
simultaneously investigating water and thermal management 
issues in a PEM fuel cell. As discussed in Ref. [40], since 
there are uncertainties concerning the GDL thermal conductiv- 
ity, parametric studies using two different thermal conductivities 
at 3 and 1.5 W(mK)~! are carried out in the present numerical 
simulations. The calculations correspond to case | in Table 3 
with the boundary temperature defined at 80°C. The temper- 
ature distributions produced using the two different thermal 
conductivities are displayed in Fig. 8. As expected, decreas- 
ing the GDL thermal conductivity increases temperature in the 
porous materials on both the anode and cathode sides. With a 
thermal conductivity at 1.5 W(m K)~!, the maximum tempera- 
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Fig. 5. Liquid saturation distributions with different liquid saturation values defined at the GDL/GC interface on the cathode side: (a) 0, (b) 0.2, and (c) 0.35. 
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Fig. 6. Effect of liquid saturation value at the GDL/GC interface on the cathode 
side on current distribution in the lateral direction. 


ture rise is more than 4°C under the gas channel, as shown in 
Fig. 8b. It should be emphasized that, in this case, the conden- 
sation and evaporation of liquid water does not contribute to the 
overall energy release but only serve to adjust the temperature 
distributions since condensation and evaporation energy com- 
pletely cancel each other, as clearly indicated in Figs. 9c and 
10b. 

The non-isothermal effect on liquid water distribution is dis- 
played in Fig. 9, in which the liquid saturation distributions 
corresponding to two different thermal conductivities are com- 
pared with that under the isothermal condition. With a thermal 
conductivity of 3W(mK)~!, the temperature increase is not 
sufficient to cause a dry region in GDL and catalyst layer on 
the cathode side although it does bring down the maximum 
liquid saturation, change liquid water distribution, and cause 
liquid water evaporation in a region near the GDL and gas chan- 
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Fig. 8. Temperature distributions under two different thermal conductivities of 
the porous materials: (a) 3 W(m K)~! and (b) 1.5 W(m K)~!. 


nel interface, as shown in Fig. 10a. The evaporation process 
occurs owing to the increased liquid saturation pressure with 
the increased temperature. As indicated in Eqs. (10) and (11), 
the evaporation rate is driven by the difference between the par- 
tial pressure of water vapor and its saturation value, and also 
controlled by the liquid saturation variation. 

As illustrated in Fig. 9c, with a thermal conductivity at 
1.5 W(mK)~!, the temperature rise in the porous materials pro- 
duces a dry region directly under the gas channel, while the water 
vapor concentration at boundary 2 remains at 15.9 molm7~?, a 
fully humidified value at 80°C. As discussed in Ref. [40], a ther- 
mal conductivity of 1.5 W(m K)™! is a typical value in carbon 
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Fig. 7. Effect of inlet humidification temperature on liquid saturation distribution at a constant cell temperature of 80°C: (a) 70°C, (b) 75°C, and (c) 80°C. 
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Fig. 9. Effect of thermal conductivity of the porous material on liquid saturation distribution: (a) isothermal case, (b) 3 W(m K)~!, and (c) 1.5 W(mK)~!. 


paper. If the thermal contact resistance is further considered, the 
liquid water distribution in Fig. 9c should represent the normal 
situations in a PEM fuel cell at the inlet region. In fact, in both the 
experimental results in Refs. [27,41], there was hardly any liq- 
uid water observed under the gas channel in the inlet region, and 
in Ref. [27], liquid droplet emergence from the GDL was only 
observable near the exit region. Further downstream along the 
cathode gas channel, with the decrease of the current density and 
consequently decreased temperature and with the water vapor in 
the gas channel becoming over-saturated, the dry region under 
the gas channel would disappear, consistent with the experimen- 
tal observations. This interesting phenomenon will be further 
investigated in the future three-dimensional numerical simula- 
tions. 

Fig. 10 shows the distributions of the condensation and evap- 
oration source terms, with a positive value for condensation and 
a negative one for evaporation. In both figures, in addition to the 
heat pipe phenomenon discussed in Ref. [20], that is evaporation 
process in the catalyst layer region (in our case the evaporation 
process is assumed to occur instantly since water is produced 
in the gaseous phase) carries out heat which is released when 
water vapor condenses near the current-collecting land, another 
heat pipe phenomenon also exists, that is evaporation process in 
the region near the gas channel carries out heat which is then 
released when water vapor condenses near the current-collecting 
land. 

Furthermore, in Fig. 10 and particularly in Fig. 10b, liquid 
water is found to be produced inside the GDL away from the cat- 
alyst layer. Since temperature is sufficiently high in the catalyst 
layer, even if water produced by the electrochemical reaction is 
assumed to be in the liquid phase, it should be evaporated and as 
aresult, only very small amount of liquid water could exist inside 
the catalyst layer at the inlet region. This is a result in favor of 
simplifying the liquid water modeling efforts in catalyst layer, 
as is the case in the present two-phase model. 


Liquid water is produced in the GDL mainly in two regions; 
one is near the land owing to the low temperature in this region 
and another further inside the GDL but still away from the cata- 
lyst layer, as shown in Fig. 10. In the open literature, the role of 
the micro-porous layer (MPL) for improving liquid flooding and 
cell performance has been discussed by a number of researchers 
based on numerical analyses of the liquid water transport pro- 
cesses [33,42]. Based on the present numerical results in Fig. 10, 
we hypothesize that at least at the inlet region of the cell, a role 
of the MPL is to prevent liquid water from entering the catalyst 
layer and therefore prevent severe liquid flooding in the region, 
since liquid water transport preferentially selects larger pores in 
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Fig. 10. Distribution of the volumetric condensation/evaporation rate under two 
different thermal conductivities of the porous materials: (a) 3 W(m K)~!, and 
(b) 1.5 W(m K)~! (unit: mol(m? s)~!). 
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the diffusion medium. Of course, this effect has not been incor- 
porated in the present two-phase PEM fuel cell model, and this 
hypothesis still needs further verification and will be considered 
in the future modeling efforts. 


4. Conclusion 


In this paper, a two-phase non-isothermal PEM fuel cell 
model, which is based on the previously developed mixed- 
domain PEM fuel cell model with a consistent treatment 
of water transport in MEA, has been established using the 
traditional two-fluid method to account for the finite-rate con- 
densation/evaporation phenomena and allow slight water vapor 
over-saturation. This two-phase PEM fuel cell model could fully 
incorporate both the anode and cathode sides, properly account 
for the various water phases, including water vapor, water in the 
membrane phase, and liquid water, and truly enable numerical 
investigations of water and thermal management issues with the 
existence of condensation/evaporation interfaces in a PEM fuel 
cell. 

This two-phase model has been applied herein in a 
two-dimensional configuration to determine the appropriate 
condensation and evaporation rate coefficients and conduct 
extensive numerical studies concerning the effects of the inlet 
humidity condition and temperature variation on liquid water 
distribution with or without a condensation/evaporation inter- 
face. The condensation rate coefficient determined in this paper, 
5000 s~!, is consistent with the detailed theoretical analyses in 
Ref. [33], and it renders the calculated water vapor concentration 
and liquid water saturation in the porous materials approach the 
thermodynamic equilibrium condition in real-world operations 
of PEM fuel cells [43]. 

Results indicate that the liquid droplet coverage model at 
the gas diffusion layer and gas channel interface exert signif- 
icant effects on liquid water distribution and cell performance 
in a PEM fuel cell. Under a low-humidity inlet condition, a 
condensation/evaporation interface would appear in the porous 
materials and its location changes with the inlet humidity value. 
Under a non-isothermal operation condition with an appropriate 
thermal conductivity of the GDL, a condensation/evaporation 
interface could also appear, resulting in a dry region in the porous 
materials directly under the gas channel, a result consistent with 
the experimental observations. 

Numerical simulations show that liquid water is mainly 
produced in the GDL in two regions; one is near the current- 
collecting land owing to the low temperature and another further 
inside the GDL but still away from the catalyst layer. This result 
suggests a new role of the micro-porous layer that is at least at 
the inlet region of the cell, it serves to prevent liquid water from 
entering the catalyst layer and thus prevent severe liquid flood- 
ing in the region. This hypothesis needs to be further verified 
and considered in the future modeling efforts. 
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